import numpy as np

class FractionalPowerPolymonial(object) :
    """
    Hyperbolic tangent kernel 
    K(x1,x2) = tanh(sc * <x1,x2> + offset) 
             =  exp( 2 * sc * <x1,x2>) - 1 / exp( ** ) + 1
    """
    def __init__(self, degree = .5, bias = 1.):
        self.degree = degree
        self.bias = bias
        
    def __call__(self, x1, x2) :
        """
        K(X, Y) = (<X, Y> + bias)^degree
        """
        K = np.dot(x1, x2)
        K += self.bias
        K = np.sign(K) *(np.abs(K) ** self.degree)
#         K_result = sqrtm(K) # K^.5 = sqrt(K)
        return K

class Kpca(object) :
    """
    Kernel PCA class 

    implemented very naively and probably unefficiently .. do not use 
    on really big dataset (>few hundreds points) 

    ex : 
    >>> k = Kpca(data) # instanciate 
    >>> k.run()        # compute eigen vectors 
    >>> p = k.project(a_vector)  # get n first projections for any vector 
    >>> p = k.projected_data()   # ....       for the input dataset

    """
    def __init__(self,dataset,kernel=GaussKernel(),components=5):
        """
        @param dataset MxN  : M row samples of dimension N 
        @param kernel is one of the kernel class (must be a callable with two vectors -> scalar ) 
        @param components : # of components keeped for projections 
        """
        self.data = dataset
        self.M = dataset.shape[0]
        self.dim = dataset.shape[1]
        self.kernel = kernel 
        self.comp = components

    def _compute_kmatrix(self):
        """ compute the Kernel dual matrix on the input dataset """
        self.K = np.empty((self.M,self.M))
        for i in range(self.M):
            self.K[i,i] = self.kernel(self.data[i],self.data[i])
            for j in range(i + 1,self.M) :
                self.K[i,j] = self.K[j,i] = self.kernel(self.data[i],self.data[j])
        print self.K
    
    def run(self) :
        """
        Compute the eigen vectors of K matrix and store them .. 
        """
        self._compute_kmatrix()
        lamb,a = np.linalg.eigh(self.K) # K is symmetric .. 
        self.alpha = []
        for i in reversed(np.argsort(lamb)[-self.comp:]) : # according to np doc , ev are not necessarily sorted .. 
            self.alpha.append(a[:,i])

    def project(self,x) :
        """
        get the projection of the vector x on the first components 
        """
        res = []
        assert x.shape == (self.dim,)
        for a in self.alpha :
            k = np.array([self.kernel(d,x) for d in self.data])
            res.append(np.dot(a,k))
        return res

    def projected_data(self) :
        """ 
        return the first components of the projected dataset 
        """
        _a = np.array(self.alpha)
        return np.dot(_a,self.K).T